#ifndef __J_PCA_H_
#define __J_PCA_H_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000

#include "DataSource.h"
#include <mkl.h>
#include <cassert>
#include <iostream>
#include <vector>

/////////////////////////////////////////////////////////////////////
//my PCA class
/////////////////////////////////////////////////////////////////////
template <typename T>
class CPCA
{
public:
	//const T *const *const m_pSamples;

	std::vector<float> m_Aver;
	std::vector<float> m_DstAver;
	std::vector<float> m_Covar;					//of size m_SrcDim*m_SrcDim
	std::vector<float> m_EigenValue;			//of size m_SrcDim
	std::vector<float> m_EigenVector;			//of size m_SrcDim*m_SrcDim

	CSimDataSource<T> *m_pResult;
private:
	int m_SrcDim;
	int m_DstDim;
	int m_nSample;
public:
	CPCA(const CDataSource<T> *const pDataSource, int dstDim=0, bool saveResult = false)
		:m_SrcDim(pDataSource->m_Dim), m_DstDim(dstDim), m_pResult(NULL)
	{
		PreProcess(pDataSource, saveResult);
	}

	CPCA(T **src, int srcDim, int nSamples, int dstDim = 0, bool saveResult = false)
		:m_SrcDim(srcDim), m_DstDim(dstDim), m_pResult(NULL)
	{
		CSimDataSource<T> dataSource(src, nSamples, srcDim);
		m_nSample = nSamples;
		PreProcess(&dataSource, saveResult);
	}

	~CPCA(){delete m_pResult;}

	void ShowData() const
	{
		std::cout.setf (std::ios::fixed | std::ios::left);
		std::cout.precision (3);

		std::cout<<"The Average Vector:\n";
		int i,j;
		for(i=0; i<m_SrcDim; i++){
			std::cout<<m_Aver[i]<<"\t";
		}

		std::cout<<std::endl<<"\nA:The CoVaraiance Matrix:\n";
		for(i=0; i<m_SrcDim; i++){
			std::cout<<i<<"Row:\t";
			for(j=0; j<m_SrcDim; j++)
				std::cout<<m_Covar[i*m_SrcDim+j]<<"\t";
			std::cout<<std::endl;
		}

		std::cout<<std::endl<<"D:The EigenValue Vector:\n";
		for(i=0; i<m_SrcDim; i++){
			std::cout<<m_EigenValue[i]<<"\t";
		}

		std::cout<<std::endl<<"\nV:The EigenVector Matrix:\n";
		for(i=0; i<m_SrcDim; i++){
			std::cout<<i<<"Row:\t";
			for(j=0; j<m_SrcDim; j++)
				std::cout<<m_EigenVector[i*m_SrcDim+j]<<"\t";
			std::cout<<std::endl;
		}

		std::vector<float> AV(m_SrcDim*m_SrcDim, 0);
		std::vector<float> VD(m_SrcDim*m_SrcDim, 0);
		std::vector<float> D(m_SrcDim*m_SrcDim, 0);
		for(i=0; i<m_SrcDim; i++){
			D[i*m_SrcDim+i] = m_EigenValue[i];
		}

		//dsymm(side, uplo, m, n, alpha, a, lda, 
		//		b, ldb, beta, c, ldc )
		const float alpha = 1.0L;
		const float beta = 0.0L;
		ssymm("L", "U", &m_SrcDim, &m_SrcDim, &alpha, &m_Covar.front(), &m_SrcDim, 
			&m_EigenVector.front(), &m_SrcDim, &beta, &AV.front(), &m_SrcDim);

		//dgemm(transa, transb, m, n, k, alpha, a, lda, 
		//		b, ldb, beta, c, ldc)
		sgemm("N", "N", &m_SrcDim, &m_SrcDim, &m_SrcDim, &alpha, &m_EigenVector.front(), &m_SrcDim, 
			&D.front(), &m_SrcDim, &beta, &VD.front(), &m_SrcDim);

		std::cout<<"\nAV:\n";
		for(i=0; i<m_SrcDim; i++){
			std::cout<<i<<"Row:\t";
			for(j=0; j<m_SrcDim; j++)
				std::cout<<AV[i*m_SrcDim+j]<<"\t";
			std::cout<<std::endl;
		}

		std::cout<<"\nVD:\n";
		for(i=0; i<m_SrcDim; i++){
			std::cout<<i<<"Row:\t";
			for(j=0; j<m_SrcDim; j++)
				std::cout<<VD[i*m_SrcDim+j]<<"\t";
			std::cout<<std::endl;
		}
	}

	void ConvertPoint(const T *const src, T *const dst)	const
	{
		int j;
		for (j=0; j<m_DstDim; j++){
			dst[j] = 0;
			for(int k=0; k<m_SrcDim; k++){
				dst[j] += (src[k] /*- m_Aver[k]*/) * m_EigenVector[k+j*m_SrcDim];  //CAUTION: when calcing the final data, do or not do the dim_mean-substraction is actually the same.
			}
		}
	}

	void ConvertData(T** src, T** dst, int nSample, T**back_data=NULL) const
	{
		for (int i=0; i<nSample; i++)
			ConvertPoint(src[i], dst[i]);
		/*show resulted PCA weight vectors*/
// 		for (int i=0; i<m_DstDim; ++i){
// 			std::cout<<"("<<i<<") : \n"; 
// 			for(int j=0; j<m_SrcDim; ++j){
// 				std::cout<<m_EigenVector[i*m_SrcDim+j]<<", "; 
// 			}
// 			std::cout<<" \n"; 
// 		}
		if(back_data==NULL) return;
		std::cout<<"\n\t get raw data back : "<<"\n"; 
		for (int i=0; i<nSample; ++i){
			for(int j=0; j<m_SrcDim; ++j){
				double sum(0.0);
				for(int k=0; k<m_DstDim; ++k){
					sum += dst[i][k] * m_EigenVector[k*m_SrcDim+j]; 
				}
				back_data[i][j] = sum+m_Aver[j];
			}
		}
	}
	void ConvertData(T** src, T* dst, int nSample, T*back_data=NULL) const
	{
		for (int i=0; i<nSample; i++)
			ConvertPoint(src[i], dst+i*m_DstDim);
	}

	void ConvertMean( )	
	{
		if(m_Aver.size()!=m_SrcDim)	return;
		m_DstAver.clear();	m_DstAver.resize(m_DstDim,0.0);
		for (int i=0; i<m_DstDim; i++){
			m_DstAver[i] = 0.0;
			for(int k=0; k<m_SrcDim; k++)
				m_DstAver[i] += m_Aver[k] * m_EigenVector[k+i*m_SrcDim];  
		}
	}

private:
	void PreProcess(const CDataSource<T> *const pDataSource, bool saveResult)
	{
		assert(pDataSource->m_Dim == m_SrcDim);
		if(pDataSource->m_Dim != m_SrcDim)
			return;

		assert(m_DstDim>0 && m_DstDim<=m_SrcDim);
		if(m_DstDim<=0 || m_DstDim>m_SrcDim)
			m_DstDim = m_SrcDim;

		const int nSamples = pDataSource->m_nSamples;
		std::cout<<"\n=====\tPCA PreProcessing with nSamples="<<nSamples<<" srcDim="<<m_SrcDim<<" dstDim="<<m_DstDim;
		clock_t tt=clock();
		//////////////////////////////////////////////////////////////////////////
		// calculating average
		//////////////////////////////////////////////////////////////////////////
		std::cout<<"\n\tCalculating Average...";
		clock_t t = clock();
		std::vector<long double> sums(m_SrcDim, 0);
		int i,j;
		std::vector<T> sample(m_SrcDim);
		for(i=0; i<nSamples; i++){
			pDataSource->getData(i, &sample.front());
			for(j=0; j<m_SrcDim; j++){
				sums[j] += sample[j];
			}
		}

		m_Aver.resize(m_SrcDim);
		for(i=0; i<m_SrcDim; i++){
			m_Aver[i] = sums[i]/nSamples;
		}

		std::vector<float> matt(nSamples*m_SrcDim);
		for(i=0; i<nSamples; i++){
			pDataSource->getData(i, &sample.front());
			for(j=0; j<m_SrcDim; j++){
				matt[i*m_SrcDim+j] = (sample[j] - m_Aver[j]);
			}
		}

		std::cout<<"\t\tfinished in "<<clock()-t<<" ms. \n";
		t=clock();

		//////////////////////////////////////////////////////////////////////////
		// calculating covariance matrix
		//////////////////////////////////////////////////////////////////////////
		std::cout<<"\n\tCalculating Covariance...";

		m_Covar.resize(m_SrcDim*m_SrcDim, 0);
		const float alpha = 1.0L;
 		const float beta = 0.0L;
		//cblas_ssyrk(CblasRowMajor, CblasUpper, CblasTrans, m_SrcDim, nSamples,
		//				1.0L, &matSample.front(), m_SrcDim, 0.0L, &m_Covar.front(), m_SrcDim);

		//ssyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc )
 		ssyrk("U", "N", &m_SrcDim, &nSamples, &alpha,
 					&matt.front(), &m_SrcDim, &beta, &m_Covar.front(), &m_SrcDim);

		matt.swap(std::vector<float>());			// cleaning
		std::cout<<"\tfinished in "<<clock()-t<<"ms";
 		t=clock();

		//////////////////////////////////////////////////////////////////////////
		// calculating eigenvalues and eigenvectors
		//////////////////////////////////////////////////////////////////////////
		std::cout<<"\n\tSolving EigenProblem...";
		m_EigenValue.resize(m_SrcDim);
		m_EigenVector.resize(m_SrcDim*m_SrcDim);
		t=clock();

		int lwork, liwork, info;
		std::vector<float> work;
		std::vector<int> iwork;

		std::vector<int> isuppz(2*m_SrcDim);
		float vl, vu, abstol=0;
		int il=m_SrcDim-m_DstDim+1, iu=m_SrcDim;
		lwork = 26*m_SrcDim;
		liwork = 10*m_SrcDim;
		work.resize(lwork);
		iwork.resize(liwork);
		matt = m_Covar;
 		//ssyevr(jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, 
		//				w, z, ldz, 
		//				isuppz, work, lwork, iwork, liwork, info)

		//ssyevr("V", "A", "U", &srcDim, &matt.front(), &m_SrcDim, &vl, &vu, &il, &iu, &abstol, &m_SrcDim, 
		ssyevr("V", "I", "U", &m_SrcDim, &matt.front(), &m_SrcDim, &vl, &vu, &il, &iu, &abstol, &m_DstDim, 
			&m_EigenValue.front(), &m_EigenVector.front(), &m_SrcDim, 
			&isuppz.front(), &work.front(), &lwork, &iwork.front(), &liwork, &info);
		
		assert(info == 0);
		if(info != 0){
			std::cout<<"Failded, Abort!"<<std::endl;
			return;
		}
		std::cout<<"resutl: "<<info<<std::endl<<" the found eigenvalues in ascending order:"<<std::endl;
//		for (unsigned i=0; i<m_EigenValue.size(); ++i)	std::cout<<m_EigenValue[i]<<", "; 
// 		//verify the orthogoanlity of between any 2 EigenVectors.
// 		for (unsigned i=0; i<m_DstDim-1; ++i){
// 			float s(.0);
// 			for(int k=0; k<m_SrcDim; ++k)
// 				s+= m_EigenVector[i+k*m_SrcDim] * m_EigenVector[i+1+k*m_SrcDim];		
// 			std::cout<<s<<" s "; 						
// 			float t(.0);
// 			for(int k=0; k<m_SrcDim; ++k)
// 				t+= m_EigenVector[i*m_SrcDim+k] * m_EigenVector[(i+1)*m_SrcDim+k];		
// 			std::cout<<t<<" t "; 						
// 		}

/*		//alternative 1: ssyevx
		lwork = 8*m_SrcDim;
		liwork = 5*m_SrcDim;
		work.resize(lwork);
		iwork.resize(liwork);
		std::vector<int> ifail(m_SrcDim);
		matt = m_Covar;
		//ssyevx(jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, 
		//		w, z, ldz
		//		work, lwork, iwork, ifail, info)
		ssyevx("V", "I", "U", &srcDim, &matt.front(), &m_SrcDim, &vl, &vu, &il, &iu, &abstol, &m_DstDim, 
			&m_EigenValue.front(), &m_EigenVector.front(), &m_SrcDim, 
			&work.front(), &lwork, &iwork.front(), &ifail.front(), &info);
		
		// alternative 2:  calculating all Eigen Values & Eigen Vectors
		j=1,i=0;
		while(j<m_SrcDim){
			j <<= 1,i++;
		}

		lwork = 3*m_SrcDim*m_SrcDim + (5 + 2*i)*m_SrcDim+1;
		liwork = 5*m_SrcDim+3;
		work.resize(lwork);
		iwork.resize(liwork);
		//dsyevd ( job , uplo , n , a , lda , w , 
		//			work , lwork , iwork , liwork , info)
		matt = m_Covar;
		ssyevd("V", "U", &m_SrcDim, &matt.front(), &m_SrcDim, &m_EigenValue.front(), 
			&work.front(), &lwork, &iwork.front(), &liwork, &info);

		//reversing;
		for(i=0; i<m_SrcDim; i++)
			std::copy(matt.begin()+i*m_SrcDim, matt.begin()+(i+1)*m_SrcDim, m_EigenVector.begin()+(m_SrcDim-i-1)*m_SrcDim);
*/

		printf("\t\tfinished in %d ms", clock()-t);

		//////////////////////////////////////////////////////////////////////////
		// All finished
		//////////////////////////////////////////////////////////////////////////
		printf("\n=====\tFinished in %d ms", clock()-tt);
		
		//ConvertMean();
		
		if(saveResult){
			printf("\n=====\tConverting source and saving result...");
			tt=clock();
			m_pResult = new CSimDataSource<T>(NULL, pDataSource->m_nSamples, m_DstDim);
			for(i=0; i<pDataSource->m_nSamples; i++){
				pDataSource->getData(i, &sample.front());
				ConvertPoint(&sample.front(), m_pResult->getData(i));
			}
			printf("finished in %d ms\n", clock()-tt);
		}
	}

};


#endif